NEW: K-Fold Validation
The first pharmacodynamic model (see fig 1) we will be looking at was produced by Friberg et. al. which models the life cycle of blood cells as transitioning through different states:
The blood cells in the proliferation state creates new cells at a rate of $k_{prol}$ when at equilibrium and transitions into the first transition state at a rate of $k_{tr}$. The blood cells also transition from one transition state to the next or from the last transition state into circulation at a rate of $k_{tr}$. The circulating blood cells are cleared from the blood stream at a rate of $k_{circ}$. However, as per the friberg paper, I will simplify the system by making $k_{prol} = k_{tr} = k_{circ} = \frac{4}{\mathrm{MTT}}$ where MTT is the mean transit time. To maintain a concentraion of $Circ$ at a constant level $Circ_0$, there is a negative feedback loop.
The drug interacts with the system by inhibiting proliferation. The relation between drug amount and drug effect is assumed to be linear, i.e. $E_\mathrm{drug} = slope*Conc_\mathrm{drug}$
This Produces equations:
$$ \dot{Prol} = k_{prol}(1-E_\mathrm{drug})(\frac{Circ_0}{Circ})^\gamma Prol - k_{tr} Prol \\ \dot{T_1} = k_{tr}Prol - k_{tr} T_1\\ \dot{T_2} = k_{tr}T_1 - k_{tr} T_2\\ \dot{T_3} = k_{tr}T_2 - k_{tr} T_3\\ \dot{Circ} = k_{tr}T_3 - k_{tr} Circ $$I initially test fit the model to simulated data before trying to fit to the real data. This shows which methods can reliably infer the parameters. The noise used had both additive and multiplicative parts, i.e. $$ X_\mathrm{Obs}(t_i) = Circ(t_i) + \left(\sigma_{base} + \sigma_{rel} Circ(t_i)^\eta\right) \epsilon_i, $$ where $$ \epsilon_i \sim \mathcal{N}\left(0,\,\sigma^{2}\right) $$ I used the typical PD parameters that Friberg et. al. inferred to create this synthesised data. The PK parameters came from the previous inference that I performed for Docetaxel. This dataset has 100 observations between the times -48 hours and 1440 hours (where dosing occurs at 0 hours).
# Options
drug = 'Simulated Drug'
dose = 2
# Retrieve the parameters for the conc-time curve
PK_params = np.load('simulated_parameters_actual_dose'+str(dose)+'.npy')
num_comp = int(len(PK_params)/2) # If this has a linear PK model
# Actual Parameters - Taken from the Friberg Paper
# (This is only to produce the simulated data. These are 'unknown' when doing the inference)
Circ_0 = 5.45
MTT = 135
gamma = 0.174
slope = 0.126
PD_actual_params = [Circ_0, MTT, gamma, slope]
param_names = ["Circ_0", "MTT", "gamma", "slope"]
# Create the Data
start_time = -48
end_time = 1440
num_obs = 100
data_times = np.linspace(start_time, end_time, num_obs)
times_before_dose = np.count_nonzero(data_times < 0)
# Conc-Time Curve
PK_times = data_times[times_before_dose:]
if PK_times[0] == 0:
conc_cuve = PK_iv_result(dose, num_comp, PK_params, PK_times)[:,0]
conc_cuve = np.concatenate((np.zeros(times_before_dose), conc_curve))
else:
conc_curve = PK_iv_result(dose, num_comp, PK_params, np.concatenate((np.zeros(1),PK_times)))[:, 0]
conc_curve = np.concatenate((np.zeros(times_before_dose), conc_curve[1:]))
# Add combined relative and constant Noise
mult_noise = np.random.normal(0, 0.453*0.1, len(data_times))
add_noise = np.random.normal(0, 0.671*0.1, len(data_times))
values_no_noise = PD_result(dose, num_comp, np.concatenate((PK_params, np.asarray(PD_actual_params))), data_times)
values_noisey = values_no_noise*(1+mult_noise)+add_noise
df = pandas.DataFrame({'TIME' : data_times, 'OBS' : values_noisey})
Before performing bayesian inference, I used optimisation to find the parameters with the highest likelihood, given the data. This should give a good starting point for using bayesian inference.
Before starting the optimisation we can estimate the parameter Circ_0.
df_before_0 = df[df["TIME"] < 0]
Circ_0_approx = sum(df_before_0["OBS"])/times_before_dose
print("approximate Circ_0: ", Circ_0_approx)
approximate Circ_0: 5.333908728293512
I also used 2 different error scores for the optimisation. Sum of squares error and the combined relative and constant gaussian noise log likelihood.
# Optimise the model with respect to the data by maximising the Log Likelihood.
problem = pints.SingleOutputProblem(pints_model_simulated, df['TIME'].to_numpy()-start_time, df['OBS'].to_numpy())
log_likelihood = ConstantAndMultiplicativeGaussianLogLikelihood(problem)
error_measure = pints.ProbabilityBasedError(log_likelihood)
lower_bound = [1e-2, 1e0, 1e-3, 1e-3, 1e-4, 1e-2, 1e-4]
upper_bound = [1e2, 1e4, 1e1, 1e1, 1e2, 1e2, 1e2]
point = np.exp((np.log(np.asarray(lower_bound)) + np.log(np.asarray(upper_bound)))/2)
point[0] = Circ_0_approx
unchanged_threshold = 1e-4
optimisation = pints.OptimisationController(error_measure, point, method=pints.CMAES, boundaries=pints.RectangularBoundaries(lower_bound, upper_bound))
optimisation.set_max_unchanged_iterations(threshold=unchanged_threshold)
optimisation.set_log_to_screen(True)
parameterslike, errorlike = optimisation.run()
actual_error = error_measure(PD_actual_params)
parameterslike = (parameterslike, errorlike, actual_error)
np.save("./Data_and_parameters/pd_sim_opt_like_params_dose_"+str(dose), parameterslike)
# Optimise the model with respect to the data
problem = pints.SingleOutputProblem(pints_model_simulated, df['TIME'].to_numpy()-start_time, df['OBS'].to_numpy())
error_measure = pints.SumOfSquaresError(problem)
lower_bound = [0.01, 0.1, 0.01, 0.001]
upper_bound = [100, 1000, 100, 10]
point = np.exp((np.log(np.asarray(lower_bound)) + np.log(np.asarray(upper_bound)))/2)
point[0] = Circ_0_approx
unchanged_threshold = 1e-4
optimisation = pints.OptimisationController(error_measure, point, method=pints.CMAES, boundaries=pints.RectangularBoundaries(lower_bound, upper_bound))
optimisation.set_max_unchanged_iterations(threshold=unchanged_threshold)
optimisation.set_log_to_screen(True)
parametersadd, erroradd = optimisation.run()
actual_error = error_measure(PD_actual_params)
parametersadd = (parametersadd, erroradd, actual_error)
np.save("./Data_and_parameters/pd_sim_opt_add_params_dose_"+str(dose), parametersadd)
# Table of Parameters
dose = 2
parameter_names =['Circ_0', 'MTT', 'gamma', 'slope', 'sigma_base', 'eta', 'sigma_rel']
parametersadd = np.load("./Data_and_parameters/pd_sim_opt_add_params_dose_"+str(dose)+".npy", allow_pickle=True)
parameterslike = np.load("./Data_and_parameters/pd_sim_opt_like_params_dose_"+str(dose)+".npy", allow_pickle=True)
PD_actual_params = np.load("./Data_and_parameters/pd_sim_actual_params_dose_"+str(dose)+".npy")
print('\t\tRun Add \t\tRun Likelihood \t\tReal')
print(parameter_names[0] + ': \t' + str(parametersadd[0][0]) + '\t' + str(parameterslike[0][0]) + '\t' + str(PD_actual_params[0]))
print(parameter_names[1] + ': \t\t' + str(parametersadd[0][1]) + '\t' + str(parameterslike[0][1]) + '\t' + str(PD_actual_params[1]))
print(parameter_names[2] + ': \t\t' + str(parametersadd[0][2]) + '\t' + str(parameterslike[0][2]) + '\t' + str(PD_actual_params[2]))
print(parameter_names[3] + ': \t\t' + str(parametersadd[0][3]) + '\t' + str(parameterslike[0][3]) + '\t' + str(PD_actual_params[3]))
print('Error: \t\t'+str(parametersadd[-1][0])+'\t'+str(parameterslike[-1][0]) )# +'\t'+str(error_measure(PD_actual_params)))
Run Add Run Likelihood Real Circ_0: 2.801241999789299 2.801241925785503 5.45 MTT: 600.7718870173026 600.7718822477289 135.0 gamma: 99.99999999999191 99.99999999999153 0.174 slope: 0.02297460272968548 0.02297460297371251 0.126 Error: 239.79480445708626 185.62452283146996
# Lets Visualise using Plotly
data_file = "./Data_and_parameters/pd_sim_data_dose_" + str(dose)
df = pandas.read_csv(data_file)
y_label = "Concentration"
x_label = "Time"
fig = px.scatter(
df,
title="Blood cell Concentration Mean",
x="TIME",
y="OBS",
width=800,
height=500,
)
fig.update_xaxes(title_text=x_label)
fig.update_yaxes(title_text=y_label)
fig.update_traces(mode='markers+lines')
fig['data'][0]['showlegend']=True
fig['data'][0]['name']='Observed Values'
fig.add_trace(go.Scatter(x=more_times, y=pints_model_simulated.simulate(parametersadd[2], more_times-start_time),
mode='lines',
name='Prediction - additive noise'))
fig.add_trace(go.Scatter(x=more_times, y=pints_model_simulated.simulate(parameterslike[2][:4], more_times-start_time),
mode='lines',
name='Prediction - combination noise'))
fig.add_trace(go.Scatter(x=more_times, y=more_values,
mode='lines',
name='Actual'))
fig.update_layout(
updatemenus=[
dict(
type = "buttons",
direction = "left",
buttons=list([
dict(
args=[{"yaxis.type": "linear"}],
label="Linear y-scale",
method="relayout"
),
dict(
args=[{"yaxis.type": "log"}],
label="Log y-scale",
method="relayout"
)
]),
pad={"r": 0, "t": -10},
showactive=True,
x=1.0,
xanchor="right",
y=1.15,
yanchor="top"
),
]
)
fig.show()
This is a useful way of estimating the correct parameters and visualising how good the model and estimated parameters fit the data. However, it would also be useful to know how sure we are on these parameter values. To do this we use Bayesian inference and MCMC to sample the parameter space and give a distribution of probable parameter values.
log_prior = pints.UniformLogPrior(lower_bound+[0]*3, upper_bound+[1000]*3)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
num_iterations = (pints_model_simulated().n_parameters()+3)*3000
startpoints = [np.concatenate((np.asarray(parametersadd), np.asarray([0.1]*3))),np.concatenate((np.asarray(parametersadd)*0.8, np.asarray([0.05]*3))), np.concatenate((np.asarray(parametersadd)*1.2, np.asarray([0.2]*3)))]
mcmc = pints.MCMCController(log_posterior, 3, startpoints, method=pints.HaarioBardenetACMC)
mcmc.set_max_iterations(num_iterations)
# mcmc.set_log_to_screen(False)
samples = mcmc.run()
pints.plot.trace(samples)
plt.show()
pints.plot.pairwise(np.vstack(samples[:,int(num_iterations/2):]))
plt.show()
The above inference was also be performed on the in vivo data provided by Roche.
import numpy as np
import math
import matplotlib.pyplot as plt
import scipy
from scipy import optimize, integrate
import pints
import pints.plot
import pandas
import plotly.express as px
import plotly.graph_objects as go
import nbformat
from ipynb.fs.full.model_simulation import PintsPDFriberg # import the pints model
from ipynb.fs.full.model_simulation import GaussianLogLikelihood
from ipynb.fs.full.model_simulation import MultiplicativeGaussianLogLikelihood
from ipynb.fs.full.model_simulation import ConstantAndMultiplicativeGaussianLogLikelihood
# Quick look at what the data contains
df = pandas.read_csv("0470-2008_2018-05-09.csv")
df = df.sort_values(['DOSE', 'TIME'], ascending=True)
group = df.groupby('DRUG')
df_view = group.apply(lambda x: x['DOSE'].unique())
df_view = df_view.apply(pandas.Series)
df_view = df_view.replace(np.nan, '', regex=True)
df_view.columns = ['Dose 1', "Dose 2", "Dose 3"]
print(df_view)
Dose 1 Dose 2 Dose 3 DRUG Controls 0.0 Docetaxel 5.0 10.0 15.0 Irinotecan 34.0 68.0 Topotecan 7.5 15.0 Vinflunine 10.0 20.0 Vinorelbine 5.0 10.0 20.0
# Options: change to one of the above drugs and corresponding dose
drug = 'Docetaxel'
observation = 'Platelets '
dose = 10
num_comp = 2
# Refine the Data
df = df.astype({'ID': 'int64'})
df_drug = df.loc[(df['DOSE'] == dose) & (df['DRUG'] == drug)] #
df_OBS = df_drug.loc[(df_drug['YNAME'] == observation)]
# Find average values
df_OBS = df_OBS.drop(df_OBS[df_OBS['OBS'] == '.'].index)[['TIME', 'DOSE', 'OBS']]
df_OBS = df_OBS.astype({'OBS': 'float64'})
df_stats = df_OBS[['TIME', 'DOSE', 'OBS']]
# df_stats = df_stats.loc[(df_stats['DOSE'] == dose)]
df_stats = df_stats.astype({'OBS': 'float64'})
df_stats = df_stats.groupby(["TIME", "DOSE"], as_index=False).filter(lambda x: len(x) > 1).groupby(["TIME", "DOSE"], as_index=False).agg({'OBS':['mean','std']})
# print(df_stats)
df_stats.columns = ['TIME', "DOSE", 'mean', 'std']
print(df_stats)
df_stats.to_csv(path_or_buf="./Data_and_parameters/pd_real_data_refined_"+drug+"_dose_"+str(dose))
# Get actual dose amount
dose_amount = df_drug.drop(df_drug[df_drug['AMT'] == '.'].index).astype({'AMT': 'float64'}).mean()['AMT']
print("Average dose amount = " + str(dose_amount))
TIME DOSE mean std 0 -48.0 10.0 931.083333 122.426570 1 24.0 10.0 773.750000 100.154464 2 48.0 10.0 676.500000 69.173694 3 72.0 10.0 711.750000 135.234056 4 96.0 10.0 642.000000 119.635558 5 120.0 10.0 808.000000 96.249675 6 144.0 10.0 780.750000 90.973897 7 168.0 10.0 711.250000 135.418795 8 192.0 10.0 756.000000 51.961524 9 216.0 10.0 997.500000 170.486559 10 240.0 10.0 1122.500000 163.752863 11 264.0 10.0 1570.333333 214.537487 12 288.0 10.0 1092.500000 320.263954 13 312.0 10.0 1443.666667 177.472064 14 336.0 10.0 1407.000000 130.761870 15 360.0 10.0 1050.333333 210.138843 16 384.0 10.0 1090.750000 75.517658 17 432.0 10.0 1047.250000 91.405233 18 480.0 10.0 826.250000 61.602895 19 552.0 10.0 946.000000 182.383662 Average dose amount = 1.9767777777777777
# OSet up the problem
problem = pints.SingleOutputProblem(pints_model_real, np.asarray(df_stats['TIME'])-start_time, df_stats['mean'])
log_likelihood = ConstantAndMultiplicativeGaussianLogLikelihood(problem, fix_eta=1)
parameter_names =['Circ_0', 'MTT', 'gamma', 'slope', 'sigma_base', 'sigma_rel']
error_measure = pints.ProbabilityBasedError(log_likelihood)
# error_measure = pints.SumOfSquaresError(problem)
lower_bound = [0.01*Circ_0_approx, 0.1, 0.01, 0.001, 0.0, 0.0]
upper_bound = [100*Circ_0_approx, 1000, 100, 10, 100, 10]
# Begin optimisation
optimisation = pints.OptimisationController(error_measure, [Circ_0_approx, 100, 1, 0.1, 0.1, 0.1], method=pints.CMAES, boundaries=pints.RectangularBoundaries(lower_bound, upper_bound))
# optimisation.set_log_to_screen(False)
parameters, error = optimisation.run()
# np.save("./Data_and_parameters/pd_real_opt_combNoEta_params_dose_"+str(dose), parameters)
# Quick visualisation
parameters = np.load("./Data_and_parameters/pd_real_opt_combNoEta_params_dose_"+str(dose)+".npy")
time_span = df_stats["TIME"].max()
times = np.linspace(start_time,time_span,1000)
# Show results
print('\t\t Parameter Estimate')
print(parameter_names[0] + ': \t' + str(parameters[0]))
print(parameter_names[1] + ': \t\t' + str(parameters[1]))
print(parameter_names[2] + ': \t\t' + str(parameters[2]))
print(parameter_names[3] + ': \t\t' + str(parameters[3]))
print('Error: \t\t'+str(error_measure(parameters)))
Run Likelihood Circ_0: 987.1821898030241 MTT: 84.29566402755697 gamma: 0.414055226474274 slope: 0.016537788114348556 Error: 143.23883610499405
# Visualisation using Plotly
y_label = drug + " Concentration"
x_label = "Time"
fig = px.scatter(
df_stats,
title=drug + " Concentration Mean",
x="TIME",
y="mean",
error_y = "std",
# facet_col="DOSE",
# color="DOSE",
width=800,
height=500,
)
fig.update_xaxes(title_text=x_label)
fig.update_yaxes(title_text=y_label)
fig.update_traces(mode='markers+lines')
fig['data'][0]['showlegend']=True
fig['data'][0]['name']='Observed Values'
fig.add_trace(go.Scatter(x=times, y=pints_model_real.simulate(parameters[:4], times-start_time),
mode='lines',
name='Prediction'))
fig.update_layout(
updatemenus=[
dict(
type = "buttons",
direction = "left",
buttons=list([
dict(
args=[{"yaxis.type": "linear"}],
label="Linear y-scale",
method="relayout"
),
dict(
args=[{"yaxis.type": "log"}],
label="Log y-scale",
method="relayout"
)
]),
pad={"r": 0, "t": -10},
showactive=True,
x=1.0,
xanchor="right",
y=1.15,
yanchor="top"
),
]
)
fig.show()
To determine the parameter distribution I performed MCMC sampling of the log posterior. However for this system, the noise term is unknown. So I will compare three different models that have the same ODE structure but different log likelihoods to capture differing noises. These noise models will be:
To compare these models I will be using WAIC and LOO-PSIS scores, which estimates the prediction accuracy for out-of-sample datapoints without requiring out-of-sample data. To do this I will need to save the pointwise loglikelihoods, i.e. $\log L(\sigma|x^{obs}_i)$ for each data point, $x^{obs}_i \in X^{obs}$, and each parameter sample, $\sigma$, obtained from the MCMC sampling.
If we assume that the distrbution of observations around the true value has both a constant and relative part, i.e. $$ X^{obs} = f\left(t, \theta\right) + \left(\sigma_{base} + \sigma_{rel} f\left(t, \theta\right)^\eta\right) \epsilon, $$ where $ \epsilon \sim N\left(0,1\right)$.
Then the likelihood of our data is given by
$$ \log L\left(\theta,\sigma_{base}, \sigma_{rel}, \eta | X^{obs}\right) = \sum_{i=1}^{n_o}\sum_{j=1}^{n_t}\left[\log L\left(\theta,\sigma_{base}, \sigma_{rel}, \eta|x^{obs}_{ij}\right)\right] $$where the pointwise log liklihoods are
$$\log L\left(\theta,\sigma, \eta|x^{obs}_{ij}\right) = - \frac{1}{2} \log 2\pi - \log \sigma_{tot,ij} -\frac{\left(x^{obs}_{ij} - f_i\left(t_j, \theta\right)\right)^2}{2\sigma_{tot,ij}^2}$$and $$\sigma_{tot,ij} = \sigma_{base, i} + f_i\left(t_j, \theta\right)^{\eta_i}\sigma_{rel, i}$$
# A look at the distribution of the parameters
parameters = np.load("./Data_and_parameters/pd_real_opt_combNoEta_params_dose_"+str(dose)+".npy")
log_likelihood = ConstantAndMultiplicativeGaussianLogLikelihoodFixEta(problem)
log_prior = pints.UniformLogPrior(np.asarray(lower_bound)*0.8, np.asarray(upper_bound)*1.2)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
num_iterations = (pints_model_real.n_parameters()+3)*3000
startpoints = [np.asarray(parameters), np.asarray(parameters)*0.8, np.asarray(parameters)*1.2]
mcmc = pints.MCMCController(log_posterior, 3, startpoints, method=pints.HaarioBardenetACMC)
mcmc.set_max_iterations(num_iterations)
# mcmc.set_log_to_screen(False)
samples = mcmc.run()
np.save("./Data_and_parameters/pd_real_mcmc_combNoEta_params_dose_"+str(dose), samples)
samples = np.load("./Data_and_parameters/pd_real_mcmc_combNoEta_params_dose_"+str(dose)+".npy")
pints.plot.trace(samples)
plt.show()
pints.plot.series(np.vstack(samples[:,int(3*num_iterations/4):]), problem)
plt.show()
pints.plot.pairwise(np.vstack(samples[:,int(3*num_iterations/4):]))
plt.show()
pints.plot.autocorrelation(np.vstack(samples[:,int(3*num_iterations/4):]))
plt.show()
If we assume that the distrbution of observations around the true value is given by a Gaussian distribution, i.e. $$ X^{obs} \sim f(t, \theta) + \sigma N(0,1), $$
then the likelihood of our data is given by
$$ \log L(\theta, \sigma | X^{obs}) = \sum_{i=1}^{n_t}\sum_{j=1}^{n_t}\left[\log L(\theta,\sigma|x^{obs}_j)\right] $$where our pointwise log liklihoods are
$$\log L(\theta,\sigma|x^{obs}_{ij}) = - \frac{1}{2} \log 2\pi - \log\sigma_i +\frac{(x^{obs}_{ij} - f_i(t_j, \theta))^2}{2\sigma_i^2}$$# A look at the distribution of the parameters
parameters = np.load("./Data_and_parameters/pd_real_opt_combNoEta_params_dose_"+str(dose)+".npy")
log_likelihood = pints.GaussianLogLikelihood(problem)
log_prior = pints.UniformLogPrior(np.asarray(lower_bound[:5])*0.8, np.asarray(upper_bound[:5])*1.2)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
num_iterations = (pints_model_real.n_parameters()+3)*3000
startpoints = [np.asarray(parameters[:5]), np.asarray(parameters[:5])*0.8, np.asarray(parameters[:5])*1.2]
mcmc = pints.MCMCController(log_posterior, 3, startpoints, method=pints.HaarioBardenetACMC)
mcmc.set_max_iterations(num_iterations)
# mcmc.set_log_to_screen(False)
samples = mcmc.run()
np.save("./Data_and_parameters/pd_real_mcmc_add_params_dose_"+str(dose), samples)
samples = np.load("./Data_and_parameters/pd_real_mcmc_add_params_dose_"+str(dose)+".npy")
pints.plot.trace(samples, parameter_names=parameter_names[:5])
plt.show()
If we assume that the distrbution of observations around the true value is relative to the magnitude of the true value, i.e. $$ X^{obs} = f\left(t, \theta\right) + \sigma f\left(t, \theta\right)^\eta \epsilon, $$ where $ \epsilon \sim N\left(0,1\right)$.
Then the likelihood of our data is given by
$$ \log L\left(\theta, \sigma, \eta | X^{obs}\right) = \sum_{i=1}^{n_t}\sum_{j=1}^{n_t}\left[\log L\left(\theta,\sigma, \eta|x^{obs}_{ij}\right)\right] $$where our pointwise log liklihoods are
$$\log L\left(\theta,\sigma, \eta|x^{obs}_{ij}\right) = - \frac{1}{2} \log 2\pi - \log f_i\left(t_j, \theta\right)^{\eta_i}\sigma_i -\frac{\left(x^{obs}_{ij} - f_i\left(t_j, \theta\right)\right)^2}{2\left(f_i\left(t_j, \theta\right)^{\eta_i}\sigma_i\right)^2}$$# A look at the distribution of the parameters
parameters = np.load("./Data_and_parameters/pd_real_opt_combNoEta_params_dose_"+str(dose)+".npy")
lower_bound = [0.01*Circ_0_approx, 0.1, 0.01, 0.001, 0.0]
upper_bound = [100*Circ_0_approx, 1000, 100, 10, 10]
log_likelihood = MultiplicativeGaussianLogLikelihoodFixEta(problem)
log_prior = pints.UniformLogPrior(np.asarray(lower_bound)*0.8, np.asarray(upper_bound)*1.2)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
num_iterations = (pints_model_real.n_parameters()+3)*3000
print(parameters)
[9.87182190e+02 8.42956640e+01 4.14055226e-01 1.65377881e-02 1.48100237e-12 1.38470015e-01]
start_point = np.concatenate((np.asarray(parameters[:4]),np.asarray(parameters[5:])))
startpoints = [start_point, start_point*0.8, start_point*1.2]
mcmc = pints.MCMCController(log_posterior, 3, startpoints, method=pints.HaarioBardenetACMC)
mcmc.set_max_iterations(num_iterations)
# mcmc.set_log_to_screen(False)
samples = mcmc.run()
np.save("./Data_and_parameters/pd_real_mcmc_multNoEta_params_dose_"+str(dose), samples)
samples = np.load("./Data_and_parameters/pd_real_mcmc_multNoEta_params_dose_"+str(dose)+'.npy')
pints.plot.trace(samples)
plt.show()
print(pints.MCMCSummary(samples))
param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------- param 1 990.76 50.73 902.33 960.69 987.82 1015.09 1133.94 1.04 1224.02 param 2 84.95 9.64 65.48 79.44 85.14 90.82 103.25 1.01 1338.61 param 3 0.43 0.12 0.22 0.36 0.42 0.49 0.71 1.00 1656.06 param 4 0.02 0.00 0.01 0.01 0.02 0.02 0.02 1.00 2550.76 param 5 0.17 0.04 0.12 0.14 0.16 0.18 0.25 1.00 901.16
import numpy as np
import pints
import pandas
# Get our Log-likelihoods:
from ipynb.fs.full.model_simulation import GaussianLogLikelihood
from ipynb.fs.full.model_simulation import MultiplicativeGaussianLogLikelihood
from ipynb.fs.full.model_simulation import ConstantAndMultiplicativeGaussianLogLikelihood
data_type = 'real'
model_type = 'multNoEta'
drug = 'Docetaxel'
observation = 'Platelets '
dose = 10
num_comp = 2
# Load data
df = pandas.read_csv("0470-2008_2018-05-09.csv")
df = df.sort_values(['DOSE', 'TIME'], ascending=True)
df = df.astype({'ID': 'int64'})
df_drug = df.loc[(df['DOSE'] == dose) & (df['DRUG'] == drug)] #
df_OBS = df_drug.loc[(df_drug['YNAME'] == observation)]
df_OBS = df_OBS.drop(df_OBS[df_OBS['OBS'] == '.'].index)[['TIME', 'DOSE', 'OBS']]
df_OBS = df_OBS.astype({'OBS': 'float64'})
df_stats = df_OBS[['TIME', 'DOSE', 'OBS']]
df_stats = df_stats.astype({'OBS': 'float64'})
df_stats = df_stats.groupby(["TIME", "DOSE"], as_index=False).filter(lambda x: len(x) > 1).groupby(["TIME", "DOSE"], as_index=False).agg({'OBS':['mean','std']})
df_stats.columns = ['TIME', "DOSE", 'mean', 'std']
# Load samples
samples = np.load("./Data_and_parameters/pd_"+data_type+"_mcmc_"+model_type+"_params_dose_"+str(dose)+".npy")
from ipynb.fs.full.model_simulation import PintsPDFriberg # import the pints model
# Create the model in PINTS
PK_params=np.load('PK_parameters_real_'+drug+str(dose)+'.npy')
start_time = df_OBS['TIME'].min()
print(start_time)
pints_model_real = PintsPDFriberg(PK_params, dose, num_comp=num_comp, start_time=start_time)
problem = pints.SingleOutputProblem(pints_model_real, np.asarray(df_stats['TIME'])-start_time, df_stats['mean'])
# Recreate
parameters = np.load("./Data_and_parameters/pd_"+data_type+"_opt_"+"combNoEta"+"_params_dose_"+str(dose)+".npy")
log_likelihood = GaussianLogLikelihood(problem)
test_likelihood = pints.GaussianLogLikelihood(problem)
num_iterations = (pints_model_real.n_parameters()+3)*3000
print(parameters)
pointwise = []
sum_likelihoods = []
no_samples = samples.shape[1]
for chain_no, chain in enumerate(samples):
proportion_check = 0
for j, sample_parameters in enumerate(chain):
single_pointwise = log_likelihood.create_pointwise_loglikelihoods(sample_parameters)
pointwise.append(single_pointwise)
test = test_likelihood.__call__(sample_parameters)
sum_likelihoods.append(test)
proportion = (j+1)/no_samples
if proportion > proportion_check:
print("chain "+str(chain_no+1)+"\t\t" +str(proportion_check*100)+"%")
if proportion_check < 0.1:
proportion_check += 0.01
else:
proportion_check += 0.1
pointwise = np.asarray(pointwise)
np.save("./Data_and_parameters/pd_"+data_type+"_mcmc_"+"add"+"_pointwiseLogLikelihoods_dose_"+str(dose)+".npy", pointwise)
Using these samples and the pointwise likelihoods I can compare the ability of these models to predict out-of-sample data. I will be using the python package Arvix to calculate the WAIC and LOO score.
import arviz as az
import matplotlib.pyplot as plt
model_type = "add"
samples = np.load("./Data_and_parameters/pd_"+data_type+"_mcmc_"+model_type+"_params_dose_"+str(dose)+".npy")
pointwise =np.load("./Data_and_parameters/pd_"+data_type+"_mcmc_"+model_type+"_pointwiseLogLikelihoods_dose_"+str(dose)+".npy")
pointwise=[pointwise[:21000,:],pointwise[21000:42000,:],pointwise[42000:,:]]
pointwise= np.asarray(pointwise)
print(pointwise.shape)
likelihood = az.convert_to_inference_data(pointwise, group='log_likelihood')
posterior = az.convert_to_inference_data(samples, group='posterior')
observed = az.convert_to_inference_data(np.asarray(df_stats['mean']), group='observed_data')
inference_data = az.concat(likelihood, posterior, observed)
I looked at the MCMC summary for each of the models, especially $\hat{r}$ to determine whether the chains converged, as well as comparing the models using both LOO-PSIS and WAIC.
model_types = ["combNoEta", "add", "multNoEta"]
compare_dict = {}
data_type = 'real'
for model_type in model_types:
samples = np.load("./Data_and_parameters/pd_"+data_type+"_mcmc_"+model_type+"_params_dose_"+str(dose)+".npy")
print(model_type)
print(pints.MCMCSummary(samples))
pointwise =np.load("./Data_and_parameters/pd_"+data_type+"_mcmc_"+model_type+"_pointwiseLogLikelihoods_dose_"+str(dose)+".npy")
pointwise=[pointwise[:21000,:],pointwise[21000:42000,:],pointwise[42000:,:]]
pointwise= np.asarray(pointwise)
likelihood = az.convert_to_inference_data(pointwise, group='log_likelihood')
posterior = az.convert_to_inference_data(samples, group='posterior')
observed = az.convert_to_inference_data(np.asarray(df_stats['mean']), group='observed_data')
inference_data = az.concat(likelihood, posterior, observed)
compare_dict[model_type] = inference_data
comparison = az.compare(compare_dict)
comparison
combNoEta param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------- param 1 983.10 48.62 873.64 957.73 985.73 1011.89 1069.35 1.03 656.58 param 2 85.26 9.19 65.07 79.64 85.43 90.75 105.03 1.01 757.39 param 3 0.44 0.12 0.24 0.36 0.43 0.50 0.75 1.01 818.39 param 4 0.02 0.00 0.01 0.01 0.02 0.02 0.02 1.01 1230.17 param 5 42.43 37.50 0.00 2.54 37.07 73.38 114.30 1.11 31.01 param 6 0.13 0.06 0.03 0.09 0.12 0.16 0.25 1.09 48.48 add param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- ------- ------- ------- ------- -------- -------- ------ ------ param 1 5652.27 5836.63 1135.26 1136.49 1363.15 14018.40 14037.24 5.26 76.46 param 2 767.39 374.51 105.92 490.60 543.92 1200.00 1200.00 3.05 102.28 param 3 6.20 8.30 0.01 0.01 1.00 18.10 18.11 4.95 57.39 param 4 2.07 4.23 0.00 0.00 0.01 1.19 12.00 4.39 112.02 param 5 0.00 0.00 0.00 0.00 0.00 0.01 0.01 5.12 98.36 multNoEta param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------- param 1 990.76 50.73 902.33 960.69 987.82 1015.09 1133.94 1.04 1224.02 param 2 84.95 9.64 65.48 79.44 85.14 90.82 103.25 1.01 1338.61 param 3 0.43 0.12 0.22 0.36 0.42 0.49 0.71 1.00 1656.06 param 4 0.02 0.00 0.01 0.01 0.02 0.02 0.02 1.00 2550.76 param 5 0.17 0.04 0.12 0.14 0.16 0.18 0.25 1.00 901.16
/home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:145: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking warnings.warn( /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:876: RuntimeWarning: overflow encountered in exp weights = 1 / np.exp(len_scale - len_scale[:, None]).sum(axis=1) /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/numpy/core/_methods.py:47: RuntimeWarning: overflow encountered in reduce return umr_sum(a, axis, dtype, out, keepdims, initial, where) /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations. warnings.warn( /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations. warnings.warn( /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:876: RuntimeWarning: overflow encountered in exp weights = 1 / np.exp(len_scale - len_scale[:, None]).sum(axis=1) /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/numpy/core/_methods.py:47: RuntimeWarning: overflow encountered in reduce return umr_sum(a, axis, dtype, out, keepdims, initial, where) /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations. warnings.warn(
| rank | loo | p_loo | d_loo | weight | se | dse | warning | loo_scale | |
|---|---|---|---|---|---|---|---|---|---|
| combNoEta | 0 | -9.506079e+01 | 4.235018e+00 | 0.000000e+00 | 1.000000e+00 | 2.765604e+00 | 0.000000e+00 | True | log |
| multNoEta | 1 | -1.318854e+02 | 4.398985e+00 | 3.682460e+01 | 0.000000e+00 | 2.864171e+00 | 2.539391e-01 | True | log |
| add | 2 | -6.836366e+29 | 6.836366e+29 | 6.836366e+29 | 6.631806e-12 | 1.657498e+29 | 1.657498e+29 | True | log |
$\hat{r}$ shows that when MCMC sampling was performed on constant gaussian noise ("add" in the table) there was no convergence in the chains. This explains the incredibly low LOO score with very high standard error. The multiplicative noise and combined noise models mostly converged and performed better in the LOO analysis with combined noise performing the best.
comparison = az.compare(compare_dict, ic="waic")
comparison
/home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:145: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking warnings.warn( /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:1405: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. See http://arxiv.org/abs/1507.04544 for details warnings.warn( /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:1405: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. See http://arxiv.org/abs/1507.04544 for details warnings.warn( /home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/arviz/stats/stats.py:1405: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. See http://arxiv.org/abs/1507.04544 for details warnings.warn(
| rank | waic | p_waic | d_waic | weight | se | dse | warning | waic_scale | |
|---|---|---|---|---|---|---|---|---|---|
| combNoEta | 0 | -9.459548e+01 | 3.769708e+00 | 0.000000e+00 | 1.000000e+00 | 2.646156e+00 | 0.000000e+00 | True | log |
| multNoEta | 1 | -1.313813e+02 | 3.894887e+00 | 3.678582e+01 | 0.000000e+00 | 2.685906e+00 | 2.268053e-01 | True | log |
| add | 2 | -2.434745e+54 | 2.434745e+54 | 2.434745e+54 | 1.800116e-12 | 1.036791e+54 | 1.036791e+54 | True | log |
WAIC also confirms the results from LOO.
For both the LOO and WAIC calculations a warning was given (LOO: Estimated shape parameter of Pareto distribution is greater than 0.7; WAIC: For one or more samples the posterior variance of the log predictive densities exceeds 0.4). This is likely due to one or more of the observations being highly influential in the fit. Thus, to confirm the results of the model comparison, I will use a more computationally intensive method, K-fold validation, that will be less impacted by these influential points.
We start this process by splitting up the data into $K$ subsets, $X^{Obs}_k$.
K = 4
model_type = "multNoEta" # choose from: "combNoEta", "add" or "multNoEta"
# Split into the K subsets
df= pandas.read_csv("./Data_and_parameters/pd_real_data_refined_"+drug+"_dose_"+str(dose))
df_shuffled = df.sample(frac=1)
data_k = np.array_split(df_shuffled, K)
Then the model needs to be fit to each set $X^{obs}_{-k} = \cup_{j \neq k} X^{obs}_j $ thus providing the posterior distribution, $P(\theta | X^{obs}_{-k})$.
# Get the data_-k for fitting
data_minus_k = [None]*K
for i in range(0, K):
data = data_k[:i] + data_k[i+1:]
data_minus_k[i] = pandas.concat(data)
data_minus_k[i] = data_minus_k[i].sort_values(['DOSE', 'TIME'], ascending=True)
# print(data_minus_k[i])
# Set up details of the MCMC run
PK_params=np.load('PK_parameters_real_'+drug+str(dose)+'.npy')
start_time = df['TIME'].min()
pints_model_real = PintsPDFriberg(PK_params, dose_amount, num_comp=num_comp, start_time=start_time)
df_before_0 = df[df["TIME"] < 0]
times_before_dose = len(df_before_0["mean"])
Circ_0_approx = sum(df_before_0["mean"])/times_before_dose
parameter_names =['Circ_0', 'MTT', 'gamma', 'slope', 'sigma_base', 'sigma_rel']
lower_bound = [0.01*Circ_0_approx, 0.1, 0.01, 0.001, 0.0, 0.0]
upper_bound = [100*Circ_0_approx, 1000, 100, 10, 100, 10]
log_posterior = [None]*K
num_iterations = (pints_model_real.n_parameters()+3)*3000
parameters = np.load("./Data_and_parameters/pd_real_opt_combNoEta_params_dose_"+str(dose)+".npy")
# Create the problem and log posteriors in Pints
for k in range(0, K):
problem = pints.SingleOutputProblem(pints_model_real, np.asarray(data_minus_k[k]['TIME'])-start_time, data_minus_k[k]['mean'])
if model_type == "add":
log_likelihood = GaussianLogLikelihood(problem)
log_prior = pints.UniformLogPrior(np.asarray(lower_bound)[:-1]*0.8, np.asarray(upper_bound)[:-1]*1.2)
startpoints = [np.asarray(parameters)[:-1], np.asarray(parameters)[:-1]*0.8, np.asarray(parameters)[:-1]*1.2]
elif model_type == "multNoEta":
log_likelihood = MultiplicativeGaussianLogLikelihood(problem, fix_eta=1)
log_prior = pints.UniformLogPrior(np.asarray(lower_bound)[:-1]*0.8, np.asarray(upper_bound)[:-1]*1.2)
startpoints = [np.asarray(parameters)[:-1], np.asarray(parameters)[:-1]*0.8, np.asarray(parameters)[:-1]*1.2]
elif model_type == "combNoEta":
log_likelihood = ConstantAndMultiplicativeGaussianLogLikelihood(problem, fix_eta=1)
log_prior = pints.UniformLogPrior(np.asarray(lower_bound)*0.8, np.asarray(upper_bound)*1.2)
startpoints = [np.asarray(parameters), np.asarray(parameters)*0.8, np.asarray(parameters)*1.2]
log_posterior[k] = pints.LogPosterior(log_likelihood, log_prior)
# run MCMC on the log posteriors
for k in range(0, K):
mcmc = pints.MCMCController(log_posterior[k], 3, startpoints, method=pints.HaarioBardenetACMC)
mcmc.set_max_iterations(num_iterations)
mcmc.set_log_to_screen(False)
samples = mcmc.run()
np.save("./Data_and_parameters/pd_real_mcmc_"+model_type+"_params_dose_"+str(dose)+"K_fold"+str(k), samples)
/home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/scipy/integrate/odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)
The out-of-sample predictive power of the model for each $x^{Obs}_i \in X^{Obs}_k$ is then given by the expected log pointwise predictive density,
$$elpd_i = \log \frac{1}{S}\sum_{s=1}^{S}P(x^{Obs}_i | \theta^{-k, s})$$where $\theta^{-k, s}$ is the $s$-th parameter sample from the MCMC sampling done with the $X^{Obs}_{-k}$ dataset and $S$ is the total number of parameters sampled. The validation score is then the sum of $elpd_i$ for all $i$.
for k in range(0, K):
# load parameter samples
samples = np.load("./Data_and_parameters/pd_real_mcmc_"+model_type+"_params_dose_"+str(dose)+"K_fold"+str(k)+".npy")
samples = np.vstack(samples[:,int(0.75*num_iterations):])
no_samples = samples.shape[0]
# create log posterior for data subset k
data_k[k] = data_k[k].sort_values(['DOSE', 'TIME'], ascending=True)
problem = pints.SingleOutputProblem(pints_model_real, np.asarray(data_k[k]['TIME'])-start_time, data_k[k]['mean'])
if model_type == "add":
log_likelihood = GaussianLogLikelihood(problem)
elif model_type == "multNoEta":
log_likelihood = MultiplicativeGaussianLogLikelihood(problem, fix_eta=1)
elif model_type == "combNoEta":
log_likelihood = ConstantAndMultiplicativeGaussianLogLikelihood(problem, fix_eta=1)
# Get pointwise likelihoods for each sample parameter
# This creates a matrix of size S x N/K
pointwise = np.zeros((no_samples, len(data_k[k])), dtype = np.longdouble)
proportion_check = 0.0
first_sample_error = 0
for j, sample_parameters in enumerate(samples):
single_pointwise = log_likelihood.create_pointwise_loglikelihoods(sample_parameters)
pointwise[j]= single_pointwise
mean = np.average(pointwise, axis=0)
pointwise_differences = np.exp(pointwise - mean) # pointwise is log likelihood, we need likelihood
elpd = - np.log(no_samples) + mean + np.log(np.sum(pointwise_differences, axis=0))
if k==0:
elpd_i = elpd
else:
elpd_i = np.concatenate((elpd_i, elpd))
# elpd_i = np.log(elpd_i)
print(elpd_i)
/home/rumney/Documents/Myleotoxicity/Myleotoxicity-PKPD/venv/lib/python3.8/site-packages/scipy/integrate/odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning) <ipython-input-34-0a2e12a2cb2f>:28: RuntimeWarning: overflow encountered in exp pointwise_differences = np.exp(pointwise - mean) # pointwise is log likelihood, we need likelihood
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf]
model_types = ["combNoEta", "add", "multNoEta"]
for model_type in model_types:
print(model_type)
for k in range(0, K):
print("k=" + str(k))
# load parameter samples
samples = np.load("./Data_and_parameters/pd_real_mcmc_"+model_type+"_params_dose_"+str(dose)+"K_fold"+str(k)+".npy")
print(pints.MCMCSummary(samples))
samples = np.vstack(samples[:,int(0.75*num_iterations):])
print("---------------------------------")
combNoEta k=0 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------- param 1 979.79 53.16 876.70 946.20 978.50 1012.34 1088.57 1.01 1064.95 param 2 81.37 12.36 57.90 72.91 81.26 89.03 106.26 1.01 917.83 param 3 0.40 0.14 0.19 0.30 0.38 0.47 0.73 1.01 991.72 param 4 0.08 0.02 0.04 0.07 0.08 0.10 0.13 1.00 1724.25 param 5 39.04 37.45 0.00 0.28 31.00 69.60 113.27 1.14 25.26 param 6 0.15 0.06 0.05 0.11 0.15 0.18 0.27 1.06 54.35 k=1 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- ------ ------ ------ ------- ------- ------- ------ ------- param 1 1003.30 59.83 865.73 971.38 1003.17 1035.07 1132.82 1.07 641.59 param 2 86.03 15.09 64.54 78.48 84.99 91.77 111.43 1.01 742.80 param 3 0.48 0.42 0.23 0.36 0.43 0.52 0.86 1.02 734.25 param 4 0.07 0.02 0.04 0.06 0.07 0.08 0.11 1.01 1521.15 param 5 35.15 36.00 0.00 0.05 24.38 62.17 112.66 1.16 26.19 param 6 0.14 0.06 0.04 0.10 0.13 0.17 0.28 1.12 46.18 k=2 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------ ------- ------ ------ param 1 841.21 143.87 663.87 725.78 774.65 968.45 1141.04 1.92 549.48 param 2 453.38 400.60 77.33 186.48 206.04 932.39 1173.07 4.00 58.79 param 3 47.04 34.17 0.36 18.13 41.93 73.66 114.44 1.18 81.36 param 4 1.57 2.97 0.03 0.04 0.05 1.37 10.39 1.56 77.30 param 5 32.69 34.01 0.00 0.85 21.69 55.96 110.38 1.16 31.92 param 6 0.18 0.11 0.04 0.10 0.14 0.25 0.42 1.86 198.69 k=3 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------- param 1 996.83 52.25 885.30 972.82 998.45 1023.38 1099.54 1.05 747.96 param 2 84.90 9.55 66.60 80.58 85.15 90.11 101.65 1.03 1078.33 param 3 0.40 0.10 0.23 0.34 0.40 0.46 0.64 1.01 1282.56 param 4 0.09 0.02 0.05 0.08 0.09 0.10 0.13 1.04 562.63 param 5 45.79 39.33 0.00 0.30 44.76 80.09 114.14 1.18 19.74 param 6 0.10 0.07 0.00 0.04 0.09 0.13 0.29 1.14 24.96 --------------------------------- add k=0 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------ param 1 723.04 409.60 51.90 539.11 791.81 1155.07 1155.07 4.47 118.51 param 2 118.98 180.51 0.12 0.12 0.31 186.01 440.05 3.89 119.10 param 3 0.04 0.04 0.01 0.01 0.05 0.06 0.06 1.21 228.13 param 4 3.80 4.55 0.01 0.01 2.15 3.92 12.00 3.10 146.48 param 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.57 118.65 k=1 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- ------ ------ ------ ------- ------- ------- ------ ------ param 1 1567.67 858.08 750.11 750.11 1348.45 2783.50 2783.50 4.66 102.26 param 2 76.53 98.73 1.40 1.40 9.99 203.93 245.32 6.77 103.22 param 3 1.32 1.45 0.10 0.10 0.54 3.17 3.48 6.72 396.88 param 4 2.76 3.98 0.00 0.00 0.03 8.59 8.59 4.97 338.45 param 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.04 120.49 k=2 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- ------ ------ ------ ------- ------- ------- ------ ------ param 1 1085.30 233.98 765.73 765.73 1176.22 1320.12 1320.12 5.86 185.72 param 2 3.26 10.47 0.29 0.29 0.61 4.60 6.90 1.05 185.91 param 3 0.06 0.07 0.01 0.01 0.02 0.12 0.12 1.51 232.28 param 4 0.59 0.49 0.01 0.01 0.55 1.21 1.21 6.79 190.19 param 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.06 184.88 k=3 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------ ------ ------ ------ ------ ------- ------- ------ ------ param 1 869.03 229.35 622.35 622.35 804.60 1174.58 1174.58 15.14 194.87 param 2 5.44 9.41 3.89 3.89 4.31 4.37 4.37 1.01 173.05 param 3 0.08 0.06 0.03 0.03 0.06 0.13 0.13 1.56 266.48 param 4 0.48 0.33 0.18 0.18 0.35 0.93 0.93 5.57 172.58 param 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.68 171.01 --------------------------------- multNoEta k=0 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- -------- -------- ------ ------ ------- -------- -------- ------ ------ param 1 15622.83 22226.69 819.55 819.55 1008.60 34893.53 68983.59 3.08 120.72 param 2 373.19 532.16 0.59 0.59 1.98 982.86 1200.00 6.19 121.57 param 3 8.40 23.38 0.01 0.01 0.02 0.06 78.82 1.55 121.33 param 4 1.22 2.61 0.00 0.00 0.02 0.02 8.19 1.60 737.51 param 5 0.00 0.01 0.00 0.00 0.00 0.01 0.02 3.11 108.27 k=1 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- ------- ------ ------ ------- -------- -------- ------ ----- param 1 6257.58 4613.27 912.45 952.34 6022.33 11554.58 14066.81 5.26 57.36 param 2 813.45 506.77 81.06 97.78 1199.58 1200.00 1200.00 4.12 65.43 param 3 6.13 14.52 0.01 0.54 3.06 10.98 13.33 1.06 62.76 param 4 3.80 5.39 0.00 0.00 0.08 10.63 12.00 6.54 30.29 param 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.40 60.94 k=2 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- -------- ------- ------- ------- -------- -------- ------ ------ param 1 8628.96 11743.03 1055.61 1209.05 1236.79 13171.78 33367.17 3.40 104.91 param 2 720.85 542.67 1.34 1.34 1047.55 1199.96 1200.00 3.43 117.52 param 3 21.04 27.58 0.01 0.29 0.38 40.84 72.22 1.82 160.10 param 4 2.61 3.40 0.00 0.04 0.04 4.59 11.21 1.90 177.83 param 5 0.01 0.02 0.00 0.00 0.00 0.00 0.05 3.92 92.75 k=3 param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ------- ------- ------- ------ ------ ------- ------- -------- ------ ----- param 1 3306.60 4251.13 853.91 998.67 1097.74 1113.80 11256.91 1.83 13.44 param 2 716.66 547.08 17.69 97.30 1127.91 1199.98 1200.00 2.22 24.94 param 3 17.90 26.91 0.01 0.05 2.12 29.38 65.85 2.57 15.71 param 4 0.72 1.33 0.00 0.01 0.05 1.19 4.67 1.44 14.34 param 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.52 10.67 ---------------------------------
We can see from $\hat{r}$ that combined noise was the only model that could converge the MCMC chains (only 1 parameter for 1 of the subsets was not close to convergence. Where as the multiplicative noise and constant noise had large difficulties in converging the chains.
# Now compare the models
# "p_waic" effective no. of parameters columns=[,"Rank", , "Weight"]
compare_K_fold = {"Model Type": ["Constant", "Relative", "Combined"]}
elpd_i = [
np.load("./Data_and_parameters/pd_real_elpd_i_K_Fold"+ "add"+"_dose_"+str(dose)+".npy"),
np.load("./Data_and_parameters/pd_real_elpd_i_K_Fold"+ "multNoEta" +"_dose_"+str(dose)+".npy"),
np.load("./Data_and_parameters/pd_real_elpd_i_K_Fold"+ "combNoEta"+"_dose_"+str(dose)+".npy")
]
compare_K_fold["ELPD"] = np.array([np.sum(elpd_i[0]), -np.sum(elpd_i[1]), np.sum(elpd_i[2])])
max_ELPD = compare_K_fold["ELPD"].max()
max_index = np.argmax(compare_K_fold["ELPD"])
compare_K_fold["d_ELPD"] = compare_K_fold["ELPD"]-min_ELPD
compare_K_fold["SE"] = np.array([np.std(elpd_i[0]), np.std(elpd_i[1]), np.std(elpd_i[2])])*np.power(len(elpd_i_comb), 0.5)
d_i = np.asarray(elpd_i) - elpd_i[max_index]
d_std = np.std(d_i, axis=1)
compare_K_fold["dSE"] = d_std*np.power(len(elpd_i_comb), 0.5)
df_compare_K_fold = pandas.DataFrame(data = compare_K_fold)
df_compare_K_fold = df_compare_K_fold.sort_values(["ELPD"], ascending=False, ignore_index=True)
df_compare_K_fold
| Model Type | ELPD | d_ELPD | SE | dSE | |
|---|---|---|---|---|---|
| 0 | Combined | -137.165095 | 0.0 | 3.276689 | 0.0 |
| 1 | Additive | -inf | -inf | NaN | NaN |
| 2 | Relative | -inf | -inf | NaN | NaN |
This lack of convergence and very small likelihoods cause numerical issues in the code meaning a K-fold comparison could not be made.